// linear-algebra.cc
//
// Program for finding dependecies from the sieving step of
// the General Number Field Sieve (GNFS).
// It does a Gaussian Elimination
// It is based on the NTL library, which is based on GnuMP

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <string>
#include <NTL/ZZ.h>
#include <NTL/RR.h>
#include <NTL/ZZX.h>
#include <NTL/ZZXFactoring.h>
#include <NTL/ZZ_pXFactoring.h>
#include <NTL/mat_GF2.h>
#include "nfsinput.h"

using namespace std;

NTL::mat_GF2* BuildMatrix(char* RFBfile, char* AFBfile, char* QCBfile, char* relations, nfsPoly poly) {
	
	NTL::ZZ a,b,norm,p,r;
	int h;
	
	int RFBnum=LineCount(RFBfile)/2; 
	std::cout << "RFB elements: " << RFBnum << std::endl; 
	int AFBnum=LineCount(AFBfile)/2;
	std::cout << "AFB elements: " << AFBnum << std::endl; 
	int QCBnum=LineCount(QCBfile)/2;
	std::cout << "QCB elements: " << QCBnum << std::endl; 
	int cols = RFBnum+AFBnum+QCBnum+1;
	//int relnum=cols+1;	
	int relnum=LineCount(relations)/2;	
	std::cout << "relations: " << relnum << std::endl; 
	std::cout << "Matrix Dimensions: " << relnum << " x " << cols << std::endl; 
	// Allocate memory for matrix
	NTL::mat_GF2* matrix = new NTL::mat_GF2(NTL::INIT_SIZE,relnum,cols);
	
	// Open files
    	std::ifstream rel(relations);
	
	// Building the matrix row by row
	for (int j=0; j<cols+1 ; j++) {
	
		// Fill in the row with 0's...
		for(int k=0;k<cols;k++) {
			matrix->put(j,k,0);
		}
		
		// Reading the relation
		rel >> a;
		rel >> b;
		norm = a+b*poly.m;
		// positive or negative norm ?
		if(norm < 0) {
			matrix->put(j,0,1);
			norm = norm * NTL::to_ZZ(-1);
		}
		
		// find and write the prime factorization of the norm
    		// ie. the contribution of the RFB
		std::ifstream rfb(RFBfile);
		h = 1;
		rfb >> p;
		rfb >> r;
		cout << "(a, b) : (" << a << ", " << b << ")" << endl;  
		while(norm != 1) {
				cout << "nb : " << norm << endl;
				cout << "(p, r) : (" << p << ", " << r << ")" << endl;
				if (norm % p == 0) {
				matrix->put(j,h,(*matrix)[j][h]+1);
				norm = norm/p;
				
			}
			else {
				rfb >> p;
				rfb >> r;
				h++;
			}
		}
    		rfb.close();
    		

		// find and write the prime factorization of the norm
    		// ie. the contribution of the AFB
		norm=NormF(poly.f,a,b);
		if (norm < 0) {
			norm = norm * -1;
		}
    		std::ifstream afb(AFBfile);
		h = 1;
		afb >> p;
		afb >> r;
		while(norm != 1) {
		                                cout << "nb : " << norm << endl;
                                cout << "(p, r) : (" << p << ", " << r << ")" << endl;

			if (norm % p == 0) {
				while ( (a+b*r)% p != 0) {
					afb >> p;
					afb >> r;
					h++;
				}
				matrix->put(j,RFBnum+h,(*matrix)[j][RFBnum+h]+1);
				norm = norm/p;
			}
			else {
				afb >> p;
				afb >> r;
				h++;
			}
		}
		afb.close();
    		
		// check for quadratric residue accordance to
		// the QCB and the Legendre symbol...
		
		
		std::ifstream qcb(QCBfile);
		for(int i=0; i < QCBnum ; i++) {
			qcb >> p;
			qcb >> r;
			norm = a + b*r;
                cout << "nb : " << norm << endl;
                                cout << "(p, r) : (" << p << ", " << r << ")" << "     ";
		cout << NTL::Jacobi(norm, p) << "    " ;			
if(NTL::Jacobi(norm,p) != 1) {
				matrix->put(j,RFBnum+AFBnum+i+1,1);
			}
		}
    		qcb.close();
	}
	// Close open files...
    	rel.close();
	cout << *matrix << endl;
	return matrix;
}


void GaussianElimination(NTL::mat_GF2* M, NTL::vec_GF2* freecols) {
	// Uses Gaussian elimination to put the matrix M
	// in reduced row-echelon form.
	// return M in reduced form and freecols will hold the 
	// free variables

	int t,h;
	bool next;
	long rows = M->NumRows();
	long cols = M->NumCols();
	NTL::vec_GF2 tmprow;
	
	h=0;
	for(int i=0;i<rows && h<cols;i++) {
		next=false;
		
		if((*M)[i][h]==0) {
			// Find row with 1 in the desired position
			t=i+1;
			while (t<rows && (*M)[t][h]==0) {
				t++;
			}
			if(t < rows) {
				// swap rows
				tmprow = (*M)[i];
				(*M)[i]=(*M)[t];
				(*M)[t]=tmprow;

			}
			else {
				// mark as free variable	
				(*freecols)[h]=1;
				i--;
				next=true;
			}
		}
		if(!next) {
			// Remove all 1's below...
			for(int j=i+1;j < rows;j++){
				if((*M)[j][h]==1) {
					(*M)[j]+=(*M)[i];
				}
			}
			// Remove all 1's above...
			for(int j=0;j < i;j++){
				if((*M)[j][h]==1) {
					(*M)[j]+=(*M)[i];
				}
			}
		}
		h++;
	}	
}	

void solveMatrix(NTL::mat_GF2* M, char* matrixfile) {

	NTL::vec_GF2 res,freecols;
	res.SetLength(M->NumRows());
	freecols.SetLength(M->NumRows());
	
	NTL::GF2 disc;
	
	NTL::mat_GF2* tmp = new NTL::mat_GF2();
	(*tmp) = NTL::transpose(*M);
	M = tmp;
	GaussianElimination(M,&freecols);
	
	std::ofstream outfile(matrixfile);
	outfile << (*M) << std::endl;
	outfile << freecols;
	outfile.close();
	std::cout << "Matrix saved in :" << matrixfile << std::endl;
}


int main (int count, char** args){
	if(count==1) {
		std::cout << std::endl;
		std::cout << "Usage: ./linear-algebra init-file"  << std::endl;
		std::cout << std::endl;
		return 1;
	}
	std::clock_t start,end;
	
	NFS_instance prob = readData(args[1]);
	std::string polyfile=args[1];
	polyfile+=".poly";
	nfsPoly poly = readPoly((char*) polyfile.c_str());

	std::string rationalfile=args[1];
	rationalfile+=".rfb";
	std::string algebraicfile=args[1];
	algebraicfile+=".afb";
	std::string quadraticfile=args[1];
	quadraticfile+=".qcb";
	std::string relationsfile=args[1];
	relationsfile+=".rel";
	std::string matrixfile=args[1];
	matrixfile+=".mat";

	NTL::RR::SetPrecision(NTL::NumBits(prob.n));
	NTL::mat_GF2* matrix;
	NTL::vec_GF2 res;
	start=std::clock();
	std::cout << "Building Matrix..." << std::endl;
	matrix=BuildMatrix((char*) rationalfile.c_str(),(char*) algebraicfile.c_str(),(char*) quadraticfile.c_str(),(char*) relationsfile.c_str(), poly);
	end=std::clock();	
	std::cout << "Matrix built...in " <<(double) (end - start)/CLOCKS_PER_SEC <<  " seconds" << std::endl;
	start=std::clock();
	std::cout << "Solving Matrix..." << std::endl;
	solveMatrix(matrix,(char*) matrixfile.c_str());
	end=std::clock();	
	std::cout << "Matrix solved...in " <<(double) (end - start)/CLOCKS_PER_SEC <<  " seconds" << std::endl;
	delete matrix;
}

